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SUMMARY 

In  the  computation  of  perturbations  for  an  orbiting  satellite,  a  problem 
can  arise  with  the  integration  of  a  rate-of-change  quantity,  the  prototype 
expression  for  which  is  k  cos  0  ,  where  0  has  a  quasi-linear  variation  with 
time.  The  problem  can  be  dealt  with  by  the  use  of  definite  integrals  between 
epochs  rather  than  indefinite  integrals  at  epochs. 
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1  INTRODUCTION  AND  BACKGROUND 

The  modelling  of  a  satellite's  orbital  motion  involves  the  integration  of 
rate-of-change  quantities  that  represent  the  effects  of  the  relevant  perturba¬ 
tions.  These  quantities  may  be  Cartesian  accelerations,  in  which  case  the 
integration  is  almost  bound  to  be  purely  numerical  and  the  general  procedure 
(method  of  'special  perturbations')1  is  usually  attributed  to  Cowell. 
Alternatively,  the  quantities  can  be  the  rates  of  change  of  a  suitable  set  of 
orbital  elements,  in  which  case  the  integration  can  still  be  entirely  numerical 
( variation-of-parameters  method)  but  is  more  likely  to  be  partly  or  wholly  based 
on  analysis  (method  of  'general  perturbations').  The  analysis  inevitably 
requires  some  approximating  assumptions,  even  when  the  rate-of-change  quantities 
have  very  simple  expressions,  the  most  important  assumption  usually  being  that 
the  different  perturbation  sources  can  be  treated  separately,  coupling  effects 
being  negligible. 

A  common  application  of  orbit  modelling  (at  RAE  in  particular)2  has  been 
the  use  of  observed  perturbations  in  the  orbits  of  one  or  more  satellites  to 
improve  the  knowledge  of  various  geophysical  parameters  associated  with  the 
atmosphere  and  the  gravitational  field3.  Typically  the  practical  analysis  is 
carried  out  in  two  stages.  In  the  first  stage  (for  each  satellite)  accumulated 
observations  of  the  satellite  itself  over  a  long  period  of  time  (perhaps  years) 
are  split  into  groups  with  each  group  spanning  a  much  shorter  period  (usually  a 
few  days) ;  each  group  of  observations  is  analysed  by  a  computer  program  (such  as 
PROP^)  that  only  has  to  model  the  dominant  perturbations,  and  in  particular  can 
ignore  variations  that  are  negligible  over  the  short  time  periods  involved.  This 
analysis  leads  to  a  series  of  element  sets,  one  set  from  each  orbit  determina¬ 
tion,  the  whole  series  spanning  the  full  period  of  the  observational  data.  In 
the  second  stage  attention  is  confined  to  a  particular  element  (eg  the 
inclination)  from  each  set,  and  its  variation,  from  set  to  set,  over  the  full 
period  is  considered.  The  object  is  to  account  for  this  variation  by  allowing 
now  for  all  the  perturbations  that  are  significant  over  the  extended  period,  and 
determining  values  (perhaps  first-time,  perhaps  revised)  for  the  numerical 
coefficients  associated  with  one  or  more  of  the  perturbation  sources  -  sometimes 
the  particular  perturbation  sources  will  have  been  omitted  from  the  first-stage 
analysis,  not  being  'dominant',  and  sometimes  it  will  be  a  matter  of  adjusting 
existing  values  that  were  adequate  for  the  first-stage  analysis. 

A  convenient  way  of  conducting  the  second-stage  analysis  is  to  remove  all 
the  known  effects  from  the  selected  orbital  element  so  that,  in  the  absence  of 
the  additional  effects  associated  with  the  perturbation  source  being  studied,  the 


Ttf  Space  384 


4 


expected  (ideal)  value  of  the  'corrected'  element  would  be  the  same  at  every 
epoch  (ie  for  each  of  the  first-stage  orbit  determinations) .  This  facilitates 
the  plotting  of  a  graph  in  which,  the  pertinent  variation  of  the  element  having 
now  been  isolated,  it  can  easily  be  judged  how  successfully  it  can  be  modelled  in 
the  evaluation  of  the  appropriate  numerical  coefficients.  A  particular  problem 
often  arises  in  this  removal,  or  'clearing',  of  known  effects,  and  it  is  this 
problem  which  is  the  subject  of  the  present  paper. 

2  THE  PROBLEM 

We  denote  by  £  the  orbital  element  assumed  to  have  been  selected  for 
analysis,  and  we  consider  the  removal  (from  the  series  of  £  values)  of  'known' 
perturbations  typified  by  the  rate-of-change  equation 


£  =  k  cos  8  . 


(1) 


The  typical  situation  here  is  as  follows:  k  is  a  slowly-changing  quantity, 
functionally  dependent  on  some  physical  parameter (s)  (such  as  a  resonant  pair  of 
harmonic  coefficients  from  the  geopotential;  the  moon's  mass;  or  a  tidal 
parameter)  and  also  on  the  three  most  constant  orbital  elements  of  the  satellite 
(a,  e  and  i);  0  ,  on  the  other  hand,  is  a  quantity  that  varies  quite  rapidly 

with  time,  typically  being  a  linear  combination  of  the  other  three  orbital 
elements  (£2,  0)  and  M)  and  quantities  such  as  sidereal  angle. 

To  identify  the  problem  associated  with  equation  (1),  we  first  suppose  that 
k  is  exactly  constant  and  that  the  rate  of  change  of  8  is  also  exactly 
constant,  equal  say  to  X  .  Then  to  clear  £  of  the  perturbations  in  question 
we  might  simply  subtract  (k/X.)  sin  0  (an  indefinite  integral)  from  the  'observed' 
value  of  £  available  at  each  epoch.  Inasmuch  as  the  values  of  k  and  X  both 
depend  on  the  values  of  the  orbital  elements,  they  could  in  principle  be  computed 
from  the  'observed'  values  of  the  element  set,  of  which  £  itself  is  one.  This 
computation  would  be  immediate  for  k  ,  but  less  so  for  X  because  it  involves  a 
rate  of  change  (eg  is  associated  with  the  even  harmonics  of  the  geopotential). 

In  practice,  neither  k  nor  X  is  constant,  but  it  is  only  the  non¬ 
constancy  of  X  that  concerns  us  so  we  will  continue  to  treat  k  as  constant. 

It  may  be  our  hope  that  the  variation  of  8  is  sufficiently  close  to  linear,  and 
the  size  of  the  perturbation  (via  the  value  of  k)  sufficiently  small  that  we  can 
escape  from  complicated  analysis,  or  more  probably  a  full-blooded  numerical 
integration,  by  just  evaluating  8  (the  local  value  of  X)  at  each  epoch  and 
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continuing  to  'correct'  each  £  by  subtracting  (k/0)  sin  0  from  it.  If  we 
trust  the  validity  of  this  approach,  the  question  then  arises  as  to  how  we  should 
calculate  0  at  each  epoch.  In  some  circumstances  (such  as  the  above-mentioned 
geopotential-defined  d>  )  a  value  can  be  computed  straightforwardly  and  with  fair 
confidence.  In  other  cases,  where  this  is  difficult,  it  may  be  thought 
preferable  to  use  a  numerical  value  for  0  ,  based  on  the  use  of  two  element 
sets.  Thus  a  value  could  be  calculated  from 

0  =  (0  -  0')/(t  -  t’)  , 

where  t'  is  the  epoch  before  (or  after)  the  t-epoch  tacitly  assumed.  The  bias 
here  is  immediately  objectionable,  and  the  use  of  epochs  both  before  and  after  t 
is  not  much  better,  being  arbitrary,  and  still  biased  unless  the  epochs  are 
evenly  spaced.  A  particularly  unsatisfactory  aspect  of  this  approach  can  be  seen 
to  arise  when  0  is  changing  very  slowly,  examples  of  this  situation  being 
associated  with  resonance  and  the  ^-induced  critical  inclination.  The  values  of 
(k/0)  sin  0  will  then  be  numerically  large  and  the  approach  eventually  loses  all 
validity . 

The  problem  may  now  be  stated  as  the  removal  of  perturbations,  of  the  type 
under  consideration,  in  such  a  way  as  to  avoid  all  difficulty  associated  with 
small  denominators,  but  without  losing  the  simplicity  attached  to  the  underlying 
assumptions.  It  should  be  noted  that  the  problem  does  not  consist  in  the  way  in 

which  a  value  might  be  estimated  for  0  ;  rather,  it  goes  back  to  the  way  in 

which  equation  (1)  should  be  integrated. 

3  THE  SOLUTION 

In  attempting  to  deal  with  equation  (1),  the  mistake  was  to  aim  at  a  series 
of  individual  quantities,  one  at  each  epoch  and  independent  of  all  the  other 
epochs,  as  the  perturbation-clearing  quantities.  This  approach,  which  led  to  the 
indefinite  integral  (k/0)  sin  0  ,  is  legitimate  for  purely  theoretical 

development,  leading  to  formulae  that  are  neat  and  compact,  but  it  is  not  at  all 

appropriate  for  the  regular  handling  of  data  in  practice. 

Let  us  suppose,  then,  that  we  have  the  N  +  1  epochs  t j  ,  where  i  =  0  , 
...,  N  ,  and  that  we  treat  0  as  constant  between  successive  epochs.  Then  the 
definite  integral  of  equation  (1)  between  epochs  t^_j  and  t^  is  given  by 
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I. 


i 


k  J  cos  8  dt 

S-i 


(k/0)  (sin  0i  -  sin  0i _ 2 )  , 


where  the  appropriate  value  of  8  is  obviously  given  by: 


0  = 


0.  -  0.  . 

i  1-1 

t.  -  t.  , 

i  1-1 


So  we  simply  write 


I . 


i 


sin  0.  -  sin  0.  , 

k  - i - ^  (t. 

0.  -  0.  ,  1 

1  1-1 


t.  J 
1-1 


(2) 


The  quantity  to  be  removed  from  £  at  t j  can  now  be  seen  to  be  I  I,  ,  if  t; 
is  regarded  as  the  'base  epoch'  at  which  no  correction  is  made.  This  may  appear 
to  introduce  a  bias  (towards  to  as  base  epoch) ,  but  the  use  of  a  different  base 
epoch  would  change  all  the  correcting  quantities,  I  I  ,  by  the  same  amount 
(independent  of  i) .  The  important  point  is  that  we  have  eliminated  zero- 
denominator  difficulties  at  a  stroke;  indeed,  equation  (2)  automatically  leads  to 

bounded  quantities,  since  I  sin  0  -  sin  cp I  <  10  -  cp|  ,  with  (in  particular) 

s 

I  I  .  <  k(tN  -  t0)  . 

It  remains  to  deal  with  a  rather  obvious  point,  namely,  that  the 
computation  of  (sin  0  -  sin  <p)  /  (0  -  <p)  has  to  be  arranged  so  as  to  avoid  rounding 
error  when  0  and  <p  are  close  together.  Since 


sin  8  -  sin  (p 

0  -  (p 


sin  S 


cos  h  (0  +  cp) 


where  8  =  h  (0  -  <p) 
to  a  polynomial  in 
criterion . 


,  this  is  a  trivial  matter  involving  the  smooth  transition3 
8  when  |8|  is  less  than  some  (possibly)  computer-dependent 
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4  FINAL  REMARKS 

The  procedure  described  is  straightforward  and,  assuming  the  avoidance  of 
rounding  error  just  referred  to,  free  from  hazard.  Some  readers  might  think  that 
the  procedure  could  be  improved  by  some  sort  of  smoothing  in  particular 
circumstances,  when  there  is  significant  'sampling  error’  in  the  observed  values 
of  0  ,  for  example.  However,  the  rewards  for  'improvements’  of  this  type  are 
likely  to  be  very  small;  at  worst,  it  is  possible  that  time  will  be  wasted  in  the 
derivation  of  corrections  that  are  actually  inferior  to  those  obtained  by  the 
unmodified  procedure. 
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